%% Compare only DMS test period data (Sensory), with Mental rotation data
%  during the response period, which should be pure mental responses.
%
% [mfr_sensory mfr_imagery]=generate_mfr_fmri_analysis( seed, movement_trial_start_times)
%
%   movement_trial_start_times (optional) tells you which trials to
%       average.  We have 5 seconds per trial.  Default: read the button
%       timing data to see if the robot moved.
function [mfr_sensory mfr_imagery]=generate_mfr_fmri_analysis2( seed, movement_trial_start_times1,...
     movement_trial_start_times2)

filename = dir(['final_seed' num2str(seed,'%03d') '_phase_4*']);
cd (filename.name);

if ~exist('movement_trial_start_times1') || isempty(movement_trial_start_times1)
    load button_timings.dat;
    movement_trial_start_times = unique(floor(button_timings(button_timings(:,2)==1,1)/5000))*5000;
else
    if size(movement_trial_start_times1,2)>1
        movement_trial_start_times1 = movement_trial_start_times1';
    end
    movement_trial_start_times = movement_trial_start_times1;
end
% Calculate mfr over movement trials only so that we can compare sensory
% and imagery conditions accurately, even with fewer movement trials during
% the sensory condition (exact matches only).
j=1;
for i = movement_trial_start_times'
    mfr_sensory_trials(:,j) = analyze_mfr_fmri(i, i+4999,.001);
    j=j+1;
end
mfr_sensory = mean(mfr_sensory_trials,2); % average over trials.

cd ..

filename = dir(['final_seed' num2str(seed,'%03d') '_phase_6*']);
cd (filename.name);

if ~exist('movement_trial_start_times2') || isempty(movement_trial_start_times2)
    load button_timings.dat;
    movement_trial_start_times = unique(floor(button_timings(button_timings(:,2)==1,1)/5000))*5000;
else
    if size(movement_trial_start_times2,2)>1
        movement_trial_start_times2 = movement_trial_start_times2';
    end
    movement_trial_start_times = movement_trial_start_times2;
end

j=1;
for i = movement_trial_start_times'
    mfr_imagery_trials(:,j) = analyze_mfr_fmri(i, i+4999,.001);
    j=j+1;
end
mfr_imagery = mean(mfr_imagery_trials,2); % average over trials.

cd ..

figure(1);
plot(mfr_sensory,mfr_imagery,'.')
axis([0 140 0 140])
axis square

%figure(2)

N_neurons = length(mfr_imagery)

active_imagery = mfr_imagery>1;
active_sensory = mfr_sensory>1;

disp('% of all neurons active in imagery');
sum(active_imagery)/N_neurons*100
disp('% of all neurons active in sensory');
sum(active_sensory)/N_neurons*100

active = active_sensory|active_imagery;
N_active = sum(active)

disp('% of ACTIVE neurons active in imagery');
sum(active_imagery)/N_active*100
disp('% of ACTIVE neurons active in sensory');
sum(active_sensory)/N_active*100


disp('% of ACTIVE active only in imagery')
sum(active_imagery&~active_sensory)/N_active*100

disp('% of ACTIVE active only in perception')
sum(active_sensory&~active_imagery)/N_active*100

disp('% overlap');
overlap = active_sensory&active_imagery; 
sum( overlap )/N_active*100

disp('corr of imagery and perception population firing rates; ACTIVE neurons')
corr(mfr_imagery(active),mfr_sensory(active))

disp('corr of imagery and perception population firing rates; overlapping neurons')
corr(mfr_imagery(overlap), mfr_sensory(overlap))

figure(2)
boxplot([mfr_imagery(overlap), mfr_sensory(overlap) ]);


disp('ranksum mfr_imagery(active) vs mfr_sensory(active):');
[p h]=ranksum(mfr_imagery(active),mfr_sensory(active))
disp('median mfr sensory:');
median(mfr_sensory(overlap))

disp('median mfr imagery:');
median(mfr_imagery(overlap))


